package ua.lviv.franko.solvers.stressanalysis.twoDimention.fem.enriched;

import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;
import java.util.HashMap;
import java.util.Iterator;
import java.util.List;
import java.util.Map;

import ua.lviv.franko.integration.Gauss2D;
import ua.lviv.franko.integration.IFunction2D;
import ua.lviv.franko.matrix.RowType;
import ua.lviv.franko.matrix.hash.lu.HashMatrixLUStyle;
import ua.lviv.franko.mlsFunctions.twoDimension.BaseFunction;
import ua.lviv.franko.mlsFunctions.twoDimension.BaseFunctionCoupled;
import ua.lviv.franko.mlsFunctions.twoDimension.Values;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.ElementMatrixForStress;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.fem.AbstractStressElement;
import ua.lviv.franko.triangulation.SaveStructures.Node;

/**
 * uses Fernandez-Mendez combining algorithm - modification of MLS shape functions.
 * 
 * @author marko.bekhta
 * 
 */
public class CombinedSolver {
	public static final int					NUMBER_OF_UNKNOWNS			= 2;
	public static final int					QUADRATURE_POINTS_NUMBER	= 7;
	public static final int					QUADRATURE_REGIONS_NUMBER	= 4;
	
	protected Gauss2D						gau;
	protected HashMatrixLUStyle					matrix;
	protected List<Double>					b;

	protected List<AbstractStressElement>	femElementList;
	protected Map<Integer, Integer>			nodeNumbersEquivalence;
	protected List<BaseFunctionCoupled>		baseFunctionsMLS;
	protected int							activeNodesNumber;

	protected double						k1;
	protected double						k2;
	protected double						k3;
	protected IFunction2D					f1Function;
	protected IFunction2D					f2Function;

	public void generate() {
		int shiftValue = activeNodesNumber * NUMBER_OF_UNKNOWNS;
		matrix = new HashMatrixLUStyle(RowType.ARRAY_TYPE, shiftValue);
		b = new ArrayList<Double>(shiftValue);
		// FEM
		for (Iterator<AbstractStressElement> iterator = femElementList.iterator(); iterator.hasNext();) {
			AbstractStressElement element = iterator.next();

			ElementMatrixForStress data = element.generateAll();
			List<Integer> pos = data.getNodeList();
			for (int i = 0; i < pos.size(); i++) {
				int firstPosition = nodeNumbersEquivalence.get(pos.get(i));
				if (firstPosition == -1) {
					continue;
				}
				for (int j = 0; j < pos.size(); j++) {
					int secondPosition = nodeNumbersEquivalence.get(pos.get(j));
					if (secondPosition == -1) {
						continue;
					}
					matrix.addValue(firstPosition * NUMBER_OF_UNKNOWNS, secondPosition * NUMBER_OF_UNKNOWNS, data.matrix[NUMBER_OF_UNKNOWNS * i][NUMBER_OF_UNKNOWNS * j]);
					matrix.addValue(firstPosition * NUMBER_OF_UNKNOWNS, secondPosition * NUMBER_OF_UNKNOWNS + 1, data.matrix[NUMBER_OF_UNKNOWNS * i][NUMBER_OF_UNKNOWNS * j + 1]);
					matrix.addValue(firstPosition * NUMBER_OF_UNKNOWNS + 1, secondPosition * NUMBER_OF_UNKNOWNS, data.matrix[NUMBER_OF_UNKNOWNS * i + 1][NUMBER_OF_UNKNOWNS * j]);
					matrix.addValue(firstPosition * NUMBER_OF_UNKNOWNS + 1, secondPosition * NUMBER_OF_UNKNOWNS + 1, data.matrix[NUMBER_OF_UNKNOWNS * i + 1][NUMBER_OF_UNKNOWNS * j
							+ 1]);
				}
				b.set(firstPosition * NUMBER_OF_UNKNOWNS, b.get(firstPosition * NUMBER_OF_UNKNOWNS) + data.vector[i * NUMBER_OF_UNKNOWNS]);
				b.set(firstPosition * NUMBER_OF_UNKNOWNS + 1, b.get(firstPosition * NUMBER_OF_UNKNOWNS + 1) + data.vector[i * NUMBER_OF_UNKNOWNS + 1]);
			}
		}

		// EFG
		for (BaseFunctionCoupled baseFunctionCoupled : baseFunctionsMLS) {
			baseFunctionCoupled.getXi().increasPosition(activeNodesNumber);
		}

		for (int i = 0; i < baseFunctionsMLS.size(); i++) {
			// System.out.println(i);
			BaseFunction baseFunction = baseFunctionsMLS.get(i);
			double[][] vals = this.gaussQuadratureIntegration(QUADRATURE_REGIONS_NUMBER, Math.max(baseFunction.getXi().getXi() - baseFunction.getXi().getRho(), -0.5),
					Math.min(baseFunction.getXi().getXi() + baseFunction.getXi().getRho(), 0.5), Math.max(baseFunction.getXi().getYi() - baseFunction.getXi().getRho(), -0.5),
					Math.min(baseFunction.getXi().getYi() + baseFunction.getXi().getRho(), 0.5), baseFunction);
			// double[][] vals = this.gaussQuadratureIntegration(-0.5, 0.5, -0.5, 0.5, baseFunction);
			int pos1 = baseFunction.getXi().getPosition();
			for (int j = 0; j < vals.length; j++) {
				int pos2 = baseFunction.getNodeList().get(j).getPosition();
				matrix.setValue(pos1 * NUMBER_OF_UNKNOWNS, pos2 * NUMBER_OF_UNKNOWNS, vals[j][0]);
				matrix.setValue(pos1 * NUMBER_OF_UNKNOWNS, pos2 * NUMBER_OF_UNKNOWNS + 1, vals[j][1]);
				matrix.setValue(pos1 * NUMBER_OF_UNKNOWNS + 1, pos2 * NUMBER_OF_UNKNOWNS, vals[j][2]);
				matrix.setValue(pos1 * NUMBER_OF_UNKNOWNS + 1, pos2 * NUMBER_OF_UNKNOWNS + 1, vals[j][3]);
			}
		}

		// for (BounderyConditionStressInterface cond : bounderyConditions) {
		// cond.applyBounderyCondition(matrix, b);
		// }
	}

	public Map<Integer, Integer> prepareNodeNumbersEquivalence(List<Node> nodes) {
		activeNodesNumber = 0;
		Map<Integer, Integer> result = new HashMap<Integer, Integer>(nodes.size());
		Collections.sort(nodes, new Comparator<Node>() {

			@Override
			public int compare(Node o1, Node o2) {
				return o1.Number - o2.Number;
			}
		});
		int number = 0;
		for (int i = 0; i < nodes.size(); i++) {
			Node node = nodes.get(i);
			if (node.IsActive) {
				result.put(node.Number, number);
				number++;
			} else {
				result.put(node.Number, -1);
			}
		}
		activeNodesNumber = number;
		return result;
	}

	protected double[][] gaussQuadratureIntegration(int n, double ax, double bx, double ay, double by, BaseFunction baseFunction) {
		double hx = (bx - ax) / n;
		double hy = (by - ay) / n;

		double[][] res = new double[baseFunction.getNodeList().size()][NUMBER_OF_UNKNOWNS * NUMBER_OF_UNKNOWNS];

		for (int i = 0; i < n; i++) {
			for (int j = 0; j < n; j++) {
				double[][] tmp = gaussQuadratureIntegration(ax + i * hx, ax + (i + 1) * hx, ay + j * hy, ay + (j + 1) * hy, baseFunction);
				for (int k = 0; k < tmp.length; k++) {
					for (int q = 0; q < tmp[0].length; q++) {
						res[k][q] += tmp[k][q];
					}
				}
			}
		}

		return res;
	}

	protected double[][] gaussQuadratureIntegration(double ax, double bx, double ay, double by, BaseFunction baseFunction) {
		// first column - baseFunction*v; second column - v*baseFunction
		double[][] result = new double[baseFunction.getNodeList().size()][NUMBER_OF_UNKNOWNS * NUMBER_OF_UNKNOWNS];
		// its a position of main node - node of base function
		int pos = baseFunction.getNodeList().indexOf(baseFunction.getXi());

		double f1 = 0.0;
		double f2 = 0.0;

		for (int j = 0; j < QUADRATURE_POINTS_NUMBER; j++) {
			double temp1 = 0.0;
			double temp2 = 0.0;
			for (int i = 0; i < QUADRATURE_POINTS_NUMBER; i++) {
				double x = (ax + bx) / 2 + (bx - ax) / 2 * gau.getX(i);
				double y = (ay + by) / 2 + (by - ay) / 2 * gau.getX(j);
				Values[] values = baseFunction.calculateAllFunctions(x, y);

				for (int k = 0; k < values.length; k++) {
					result[k][0] += gau.getW(j) * gau.getW(i) * (values[pos].valX * values[k].valX * k1 + values[pos].valY * values[k].valY * k3);
					result[k][1] += gau.getW(j) * gau.getW(i) * (values[pos].valY * values[k].valX * k2 + values[pos].valX * values[k].valY * k3);
					result[k][2] += gau.getW(j) * gau.getW(i) * (values[pos].valY * values[k].valX * k2 + values[pos].valX * values[k].valY * k3);
					result[k][3] += gau.getW(j) * gau.getW(i) * (values[pos].valX * values[k].valX * k3 + values[pos].valY * values[k].valY * k1);

				}

				temp1 += gau.getW(i) * values[pos].val * f1Function.calculate(x, y);
				temp2 += gau.getW(i) * values[pos].val * f2Function.calculate(x, y);
			}
			f1 += temp1 * gau.getW(j);
			f2 += temp2 * gau.getW(j);
		}
		double coef = (bx - ax) * (by - ay) / 4.0;
		f1 *= coef;
		f2 *= coef;
		for (int i = 0; i < result.length; i++) {
			result[i][0] *= coef;
			result[i][1] *= coef;
			result[i][2] *= coef;
			result[i][3] *= coef;
		}
		b.set(baseFunction.getXi().getPosition() * NUMBER_OF_UNKNOWNS, b.get(baseFunction.getXi().getPosition() * NUMBER_OF_UNKNOWNS) + f1);
		b.set(baseFunction.getXi().getPosition() * NUMBER_OF_UNKNOWNS + 1, b.get(baseFunction.getXi().getPosition() * NUMBER_OF_UNKNOWNS + 1) + f2);
		return result;
	}
}
